
10/500369 



DESCRIPTION 

Brain Current Source Estimating Method, Brain Current Source Estimating Program, 
Recording Medium Recording Brain Current Source Estimating Program and Brain 
5 Current Source Estimating Apparatus 

Technical Field 

The present invention relates to a method of estimating position or distribution 
of a brain current source or sources generating magnetic fields or electric fields outside a 
10 scalp in correspondence to brain activities, as well as to a brain current source 

estimating program, a recording medium recording the brain current source estimating 
program and to a brain current source estimating apparatus. 



Background Art 

15 Thanks to remarkable development in the field of biomedical measurement 

technique, precision in measuring weak electric field (brain wave) or weak magnetic 
field (brain magnetic wave) generated from the brain, of which measurement has been 
difficult and error-prone, has been improved year by year. 

Specifically, receiving external stimuli, neural cells in the brain generate a current. 

20 The current results in the afore-mentioned weak electric field or weak magnetic field. 
Here, "brain wave" refers to the electric field generated from the brain by the current 
from the neural cells. An "electroencephalogram: EEG" refers to a method of 
measuring the brain wave. 

The "brain magnetic wave" refers to a magnetic field generated from the brain by 

25 the current from the neural cells. A "magnetoencephalography: MEG" refers to a 
method of measuring the brain magnetic wave. A crucial advantage of the 
magnetoencephalography is that, as the magnetic field is almost free of any influence of 
volume conductor, it is expected that relatively accurate three-dimensional estimation of 
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the position of a brain current source can be attained by measuring magnetism from 
outside one's scalp. 

In the analysis of the brain magnetic wave, an active portion of the brain is 
estimated in a non-invasive manner, by measuring the generated magnetic field from 
5 outside the brain. 

The magnetic field, however, is so weak that it is very susceptible to the 
influence of external magnetic field such as terrestrial magnetism. Therefore, the weak 
magnetic field is measured by a Superconducting QUantum Interference Device 
(SQUID), which is a measuring device utilizing superconductivity, within a shield that 
10 shuts out any external magnetic field. 

It is noted, however, that in the field of studying algorithms for "estimating 
positions of brain current source," decisive method is non-existent at present, though 
various and many variations of initial models have been tried. 

By way of example, "dipole estimation method" as one algorithm for "estimating 
15 positions of brain current sources," is disclosed in Reference 1 : J. C. Mosher, P. S. 

Lewis and R. M. Leahy, IEEE Trans. Biomed. Engng. <1992> vol. 39, pp.541-557. In 
the "dipole estimation method," however, the position of a dipole is estimated from 
observed magnetic field, assuming that the current source in the brain can be represented 
by one or a number of current dipoles, and this method is disadvantageous in that it is 
20 difficult to determine the number of dipoles. 

As another algorithm, "spatial filtering method" is disclosed in Reference 2: K. 
Toyama, K. Yoshikawa, Y. Yoshida, Y. Kondo, S. Tomita, Y. Takanashi, Y. Ejima and 
S. Yoshizawa, Neuroscience, 1999, 91 (2), pp. 405-415. In the "spatial filtering 
method," location of a brain current source is restricted in consideration of physiological 
25 findings, and distribution of dipoles is estimated. This method is disadvantageous in 
that accurate estimation of the depth of the current source is not possible. 

Disclosure of the Invention 
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An object of the present invention is to provide a brain current source estimating 
method that enables estimation of a position, depth direction inclusive, of a brain current 
source. 

Another object of the present invention is to provide a brain current source 
5 estimating method that enables estimation of positions of a plurality of brain current 
sources from observed data, even when there are a plurality of brain current sources. 

A still further object of the present invention is to provide a brain current source 
estimating method enabling estimation with higher accuracy, when observation data 
obtained by a method of observation independent of the observation of magnetic field 
10 for estimating the brain current source are available, by combining data obtained by the 
plurality of observation methods. 

A still further object of the present invention is to provide a brain current source 
estimating program that enables estimation of a position, depth direction inclusive, of a 
brain current source and a recording medium having the program recorded thereon. 
15 A still further object of the present invention is to provide a brain current source 

estimating program that enables estimation of positions of a plurality of brain current 
sources from observed data, even when there are a plurality of brain current sources, 
and a recording medium having the program recorded thereon. 

A still further object of the present invention is to provide a brain current source 
20 estimating program enabling estimation with higher accuracy, when observation data 
obtained by a method of observation independent of the observation of magnetic field 
for estimating the brain current source are available, by combining data obtained by the 
plurality of observation methods, and a recording medium having the program recorded 
thereon. 

25 A still further object of the present invention is to provide a brain current source 

estimating apparatus that enables estimation of a position, depth direction inclusive, of a 
brain current source. 

A still further object of the present invention is to provide a brain current source 
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estimating apparatus that enables estimation of positions of a plurality of brain current 
sources from observed data, even when there are a plurality of brain current sources. 

A still further object of the present invention is to provide a brain current source 
estimating apparatus enabling estimation with higher accuracy, when observation data 
5 obtained by a method of observation independent of the observation of magnetic field 
for estimating the brain current source are available, by combining data obtained by the 
plurality of observation methods. 

In order to attain the above described objects, the present invention provides a 
brain current source estimating method for estimating, based on an electromagnetic field 

10 observed outside a scalp, a position of a current source as a source of the 

electromagnetic wave existing in the brain, including the steps of: setting, in the brain, a 
plurality of virtual curved surfaces having depths from brain surface different from each 
other and shapes not intersecting with each other and setting lattice points on each of 
the virtual curved surfaces; estimating, on each of the virtual curved surfaces, a current 

15 distribution for recovering the observed electromagnetic field; and based on an 

expansion of the current distribution estimated on the virtual curved surfaces and a 
difference between the electromagnetic field recovered based on the current distribution 
and the observed electromagnetic field, identifying a virtual curved surface at which the 
expansion and the difference attain relative minimums among the plurality of virtual 

20 curved surfaces as a true curved surface on which the current source exists. 

Preferably, the step of estimating the current distribution includes the step of 
determining posterior probability by Bayesian estimation method from prior distribution 
and observation data of the electromagnetic field; and the step of identifying as a true 
curved surface on which the current source exists includes the step of identifying a 

25 virtual curved surface of which corresponding posterior probability attains the highest, 
among the virtual curved surfaces. 

Preferably, the step of estimating a current distribution includes the step of 
identifying a first virtual curved surface closest to the brain surface and having posterior 
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probability attaining a relative maximum, among the plurality of virtual surfaces, while 
successively moving from a virtual curved surface on the side of the brain surface to a 
deeper side; and the step of identifying a curved surface as a true curved surface on 
which the current source exists includes the steps of identifying a localized current 
5 distribution corresponding to a point of relative maximum of the current distribution, on 
the first virtual curved surface, separating a plurality of local surfaces each including the 
localized current distribution, and fixing, among the plurality of local surfaces, local 
surfaces other than a local surface as an object of identification, moving the local surface 
as an object of identification in the depth direction, and identifying positions where the 
10 posterior probability attains the relative maximum, successively from the side closer to 
the brain surface. 

Preferably, in the step of estimating a current distribution, the current 
distribution is estimated with a first spatial resolution; and the method further includes 
the step of re-estimating the current distribution with a second spatial resolution higher 

15 than the first resolution and resolution of the plurality of virtual curved surfaces in the 
depth direction being improved. 

Preferably, the step of estimating a current distribution includes the step of 
setting, when the current distribution is estimated in accordance with Bayesian 
estimation, a hierarchical prior distribution in the Bayesian estimation using observation 

20 data obtained by other observation method independent of the observation of 

electromagnetic field for the estimation of the current source. By way of example, 
when it is known from observation data obtained by a different method of observation 
that the location of the brain current source is within a restricted area, search in the 
depth direction may be omitted and the current distribution in the restricted area can be 

25 obtained by Bayesian estimation. 

According to another aspect, the present invention provides a program for a 
computer for estimating, based on an electromagnetic field observed outside a scalp, a 
position of a current source as a source of the electromagnetic wave existing in the brain, 
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to have the computer execute the steps of: setting, in the brain, a plurality of virtual 
curved surfaces having depths from brain surface different from each other and shapes 
not intersecting with each other and setting lattice points on each of the virtual curved 
surfaces; estimating, on each of the virtual curved surfaces, a current distribution for 
5 recovering the observed electromagnetic field; and based on an expansion of the current 
distribution estimated on the virtual curved surfaces and a difference between the 
electromagnetic field recovered based on the current distribution and the observed 
electromagnetic field, identifying a virtual curved surface at which the expansion and 
the difference attain relative minimums among the plurality of virtual curved surfaces 

10 as a true curved surface on which the current source exists. 

Preferably, the step of estimating the current distribution includes the step of 
determining posterior probability by Bayesian estimation method from prior distribution 
and observation data of the electromagnetic field; and the step of identifying as a true 
curved surface on which the current source exists includes the step of identifying a 

15 virtual curved surface of which corresponding posterior probability attains the highest, 
among the virtual curved surfaces. 

Preferably, the step of estimating a current distribution includes the step of 
identifying a first virtual curved surface closest to the brain surface and having posterior 
probability attaining a relative maximum, among the plurality of virtual surfaces, while 

20 successively moving from a virtual curved surface on the side of the brain surface to a 
deeper side; and the step of identifying a curved surface as a true curved surface on 
which the current source exists includes the steps of identifying a localized current 
distribution corresponding to a point of relative maximum of the current distribution, on 
the first virtual curved surface, separating a plurality of local surfaces each including the 

25 localized current distribution, and fixing, among the plurality of local surfaces, local 

surfaces other than a local surface as an object of identification, moving the local surface 
as an object of identification in the depth direction, and identifying positions where the 
posterior probability attains the relative maximum, successively from the side closer to 
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the brain surface. 

Preferably, in the step of estimating a current distribution, the current 
distribution is estimated with a first spatial resolution; and the method further includes 
the step of re-estimating the current distribution with a second spatial resolution higher 
than the first resolution and resolution of the plurality of virtual curved surfaces in the 
depth direction being improved. 

Preferably, the step of estimating a current distribution includes the step of 
setting, when the current distribution is estimated in accordance with Bayesian 
estimation, a hierarchical prior distribution in the Bayesian estimation using observation 
data obtained by other observation method independent of the observation of 
electromagnetic field for the estimation of the current source. By way of example, 
when it is known from observation data obtained by a different method of observation 
that the location of the brain current source is within a restricted area, search in the 
depth direction may be omitted and the current distribution in the restricted area can be 
obtained by Bayesian estimation. 

According to a still further aspect, the present invention provides a brain current 
source estimating apparatus for estimating, based on an electromagnetic field observed 
outside a scalp, a position of a current source as a source of the electromagnetic wave 
existing in the brain, including: virtual curved surface setting means for setting, in the 
brain, a plurality of virtual curved surfaces having depths from brain surface different 
from each other and shapes not intersecting with each other and setting lattice points on 
each of the virtual curved surfaces; current distribution estimating means for estimating, 
on each of the virtual curved surfaces, a current distribution for recovering the observed 
electromagnetic field; and current source identifying means for identifying, based on an 
expansion of the current distribution estimated on the virtual curved surfaces and a 
difference between the electromagnetic field recovered based on the current distribution 
and the observed electromagnetic field, a virtual curved surface at which the expansion 
and the difference attain relative minimums among the plurality of virtual curved 



surfaces as a true curved surface on which the current source exists. 

Preferably, the current distribution estimating means includes posterior 
probability determining means for determining posterior probability by Bayesian 
estimation method from prior distribution and observation data of the electromagnetic 
5 field; and the current source identifying means includes virtual curved surface identifying 
means for identifying a virtual curved surface of which corresponding posterior 
probability attains the highest, among the virtual curved surfaces. 

Preferably, the current distribution estimating means includes shallowest virtual 
curved surface identifying means for identifying a first virtual curved surface closest to 

10 the brain surface and having posterior probability attaining a relative maximum, among 
the plurality of virtual surfaces, while successively moving from a virtual curved surface 
on the side of the brain surface to a deeper side; and the current source identifying 
means includes localized current distribution identifying means for identifying a localized 
current distribution corresponding to a point of relative maximum of the current 

15 distribution, on the first virtual curved surface, local surface extracting means for 

separating a plurality of local surfaces each including the localized current distribution, 
and local surface position identifying means for fixing, among the plurality of local 
surfaces, local surfaces other than a local surface as an object of identification, moving 
the local surface as an object of identification in the depth direction, and identifying 

20 positions where the posterior probability attains the relative maximum, successively from 
the side closer to the brain surface. 

Preferably, the current distribution estimating means estimates the current 
distribution with a first spatial resolution and thereafter re-estimates the current 
distribution with a second spatial resolution higher than the first resolution and 

25 resolution of the plurality of virtual curved surfaces in the depth direction being 
improved. 

Preferably, the current distribution estimating means includes means for setting, 
when the current distribution is estimated in accordance with Bayesian estimation, a 



-8- 



hierarchical prior distribution in the Bayesian estimation using observation data obtained 
by other observation method independent of the observation of electromagnetic field for 
the estimation of the current source. By way of example, when it is known from 
observation data obtained by a different method of observation that the location of the 
5 brain current source is within a restricted area, search in the depth direction may be 

omitted and the current distribution in the restricted area can be obtained by Bayesian 
estimation. 

According to the brain current source estimating method of the present invention, 
it is possible to estimate the position of a brain current source including the depth 
10 direction, from observation data of MEG or EEG. Such estimation in the depth 

direction is applicable even when there are a plurality of current sources. Further, the 
method is applicable where the current sources exist locally as in the case of current 
dipoles and applicable to a current source that has a certain expansion. Further, it is 
possible to estimate how the current source expands. 

15 

Brief Description of the Drawings 

Fig. 1 is a schematic illustration of an exemplary configuration of an MEG 

system. 

Fig. 2 is a schematic illustration representing a magnetic field generated by a 
20 current source, observed on an appropriate curved surface. 

Fig. 3 is a schematic illustration representing a relation between a current source 
in the brain and a plurality of virtual curved surfaces. 

Fig. 4 is a flow chart representing an overall procedure of the brain current 
source estimating method in accordance with the present invention. 
25 Fig. 5 is a flow chart representing a process of initial estimation of the current 

source using an initial resolution. 

Fig. 6 is a flow chart representing a process for identifying a current source 
closest to the brain surface. 
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Fig. 7 is a first flow chart representing a process for identifying depth of a 
current source corresponding to each local surface. 

Fig. 8 is a second flow chart representing a process for identifying depth of a 
current source corresponding to each local surface. 
5 Fig. 9 is a first flow chart representing a process for re-estimating a position of a 

current source with higher spatial resolution. 

Fig. 10 is a second flow chart representing a process for re-estimating a position 
of a current source with higher spatial resolution. 

Fig. 1 1 is a top view of a magnetic field distribution observed on a surface of a 
10 hemisphere, assuming that the human brain is a hemisphere. 

Figs. 12(a), 12(b) and 12(c) represent results of initial estimation of current 
sources using initial resolution, where the radius is R=5.0, R=6.0 and R=7.0, 
respectively. 

Figs. 13(d), 13(e) and 13(f) represent results of initial estimation of current 
15 sources using initial resolution, where the radius is R=7.0, R=8.0 and R=9.0, 
respectively, that is, closer to the surface of the brain. 

Fig. 14 represents magnitude of free energy on each virtual hemisphere, obtained 
when current distribution that attains highest free energy is calculated. 

Figs. 15(a), 15(b) and 15(c) represent processes for identifying the current 
20 source executed further on a local surface including a point of relative maximum, where 
the radius is R=5.0, R=6.0 and R=7.0, respectively. 

Figs. 16(d), 16(e) and 16(f) represent processes for identifying the current 
source executed further on a local surface including a point of relative maximum, where 
the radius is R=7.5, R=8.0 and R=9.0, respectively. 
25 Fig. 17 represents the free energy calculated with the depth of local surface 

varied. 

Fig. 18 is a top view of a magnetic field distribution observed on a surface of the 
human brain assumed as a hemisphere, when there are two current sources in the brain. 
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Figs. 19(a), 19(b) and 19(c) represent results of initial estimation of current 
sources using initial resolution, where the radius is R=5.0, R=6.0 and R=7.0, 
respectively. 

Figs. 20(d), 20(e) and 20(f) represent results of initial estimation of current 
5 sources using initial resolution, where the radius is R=7.5, R=8.0 and R=9.0, 
respectively. 

Fig. 21 represents radius dependency of free energy, when current distribution 
that attains highest free energy is calculated for each virtual hemispherical surface. 

Figs. 22(a), 22(b) and 22(c) represent current distribution on local surfaces 
10 where the radius is R=5.0, R=6.0 and R=7.0, respectively, illustrating the process for 
calculating the depth of a first local surface attaining maximum posterior probability, to 
identify a current source closest to the brain surface. 

Figs. 23(d), 23(e) and 23(f) represent current distribution on local surfaces 
where the radius is R=7.5, R=8.0 and R=9.0, respectively, illustrating the process for 
15 calculating the depth of a first local surface attaining maximum posterior probability, to 
identify a current source closest to the brain surface. 

Fig. 24 represents local surface position (radius) dependency of free energy 
when a virtual local surface is moved. 

Figs. 25(a), 25(b) and 25(c) represent current distribution on local surfaces 
20 where the radius is R=5.0, R=5.5 and R=6.0, when one local surface is fixed on a 

specific surface and a local surface corresponding to the other current source is moved. 

Figs. 26(d), 26(e) and 26(f) represent current distribution on local surfaces 
where the radius is R=6.5, R=7.0 and R=7.5, when. one local surface is fixed on a 
specific surface and a local surface corresponding to the other current source is moved. 
25 Fig. 27 represents local surface position (radius) dependency of free energy 

when a virtual local surface is moved. 

Fig. 28 is a first graph representing a result of simulation considering localized 
condition only. 
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Fig. 29 is a second graph representing a result of simulation considering 
localized condition only. 

Fig. 30 is a first graph representing a result of simulation considering both 
continuous condition and localized condition. 
5 Fig. 3 1 is a second graph representing a result of simulation considering both 

continuous condition and localized condition. 

Fig. 32 is a first graph representing a result of simulation considering both 
continuous condition and localized condition, and using hierarchical prior distribution 
that facilitates estimation of a current in a region where presence of a current source is 
10 highly likely. 

Fig. 33 is a second graph representing a result of simulation considering both 
continuous condition and localized condition, and using hierarchical prior distribution 
that facilitates estimation of a current in a region where presence of a current source is 
highly likely. 

15 

Best Modes for Carrying Out the Invention 

Embodiments of the present invention will be described with reference to the 

figures. 

As already described, magnetoencephaography: MEG and 
20 electroencephalogram: EEG have been known as methods of monitoring brain activities 
from the outside. In the following, description will be given mainly focusing on MEG. 
It is noted, however, that the present invention is also applicable to EEG, when the 
magnetic field in MEG is replaced by the electric field. 

The term "electromagnetic field" originally refers to co-existence of "electric 
25 field" and "magnetic field." In the present specification, however, an expression 

"observe an electromagnetic field" will be used generally, where the physical amount "to 
be observed" is an "electric field" or a "magnetic field." 

[First embodiment] 
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Fig. 1 is a schematic illustration of an exemplary configuration of an MEG 

system. 

An MEG 12 includes an array of multi-channel SQUID flux meter (super 
sensitive magnetometer), for measuring a magnetic field on the scalp surface of a subject 
5 10. Receiving the result of measurement by MEG 12, a computer system 20 analyses 
the result of measurement, and performs a process for estimating the position of a brain 
current source. 

The present invention relates to software processing by computer system 20. It 
is noted, however, that part of or all of the process may be performed by hardware to 

10 increase the speed of processing. The software is not specifically limited, and it may be 
recorded on a recording medium 22 such as a CD-ROM (Compact Disk Read Only 
Memory) or a DVD-ROM (Digital Versatile Disc Read Only Memory) and installed in 
computer system 20, or it may be distributed through a communication network. 
[Principle of brain current source estimation] 

15 Prior to detailed description of the "brain current source estimating method" in 

accordance with the present invention, the principle and outline of the estimation 
method of the present invention will be summarized. 

(Recovery of electromagnetic field by current distribution on a curved surface) 
It is well-known as the principle of electromagnetic shield that when a current 

20 source is surrounded by an ideal electromagnetic shielding surface formed of metal and 
magnetic body, electromagnetic field cease to exist outside the shielding surface. 
Specifically, by the electromagnetic field formed by the current flowing over the 
shielding surface and a collection of small magnets constituting the magnetic body, the 
electromagnetic field formed by the current source outside the shielding surface is fully 

25 cancelled. Further, the small magnets can be replaced by a virtual eddy current. 

When viewed from a different point, this means that an electromagnetic field identical to 
that formed by the current source can be formed outside the shielding surface by causing 
an appropriate current flow over the shielding surface. On the contrary, the electric 
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field outside the shielding surface cannot be fully recovered, whatever current is caused 
to flow over the shielding surface. This fact will be referred to as the "principle of 
electromagnetic field recovery." 

Fig. 2 is a schematic illustration representing a magnetic field generated by a 
5 current source, observed on an appropriate curved surface. 

As can be seen from Fig. 2, when a virtual surface is prepared between an 
observing surface and the current source, it is possible to recover the electromagnetic 
field formed by the current source on the observing surface by causing an appropriate 
flow of current on the virtual surface, in accordance with the "principle of 
10 electromagnetic field recovery." 

Fig. 3 is a schematic illustration representing a relation between a current source 
in the brain and a plurality of virtual curved surfaces. 

Referring to Fig. 3, considering that the electromagnetic field formed by the 
current attenuates in inverse proportion to the square of distance, the expansion of 
15 current on the virtual curved surface (virtual curved surface 2) equivalent to the current 
source becomes wider as the virtual surface is further away from the current source. 
Therefore, the expansion of the current on virtual curved surface 1 is wider than that on 
virtual curved surface 2. 

On a virtual curved surface (virtual curved surface 3) on the side opposite to the 
20 observing surface with respect to the current source, it is impossible to fully recover the 
electromagnetic field formed on the observing surface by the current source. 

According to the present invention, based on the principle described above, the 
current source is estimated from the observed data of the electromagnetic field on the 
observing surface. 
25 (Principle of current source estimation) 

Assume that a magnetic field (or an electric field) formed by a current source 
generated in the brain is observed on an observing surface that is close to the surface of 
the brain, as shown in Fig. 2. 



- 14- 



A virtual curved surface in the brain is considered, and current distribution on 
the virtual curved surface that recovers the observed magnetic field is calculated. 
When the virtual curved surface is moved to the inner side of the brain with the radius 
made gradually smaller, the expansion of the current distribution recovering the 
5 observed magnetic becomes smaller, and it becomes the smallest when the virtual 
surface encompasses the true current source. When the virtual curved surface is 
further moved to be deeper than the current source, the expansion of current distribution 
comes to be wider again, and the difference between the magnetic field generated by the 
current and the observed magnetic field also becomes larger. 

10 Accordingly, the depth of the current source can be identified by reviewing the 

expansion of the current distribution recovering the observed magnetic field and the 
error in recovery of the magnetic field. Further, by calculating the current distribution 
on the virtual curved surface at the depth identified in this manner, the expansion of the 
current source can also be found. The foregoing is the principle of the present 

15 invention. 

(Identification of current source depth by Bayesian estimation) 
In order to specifically implement the principle of current source estimation 
described above, in the present invention, a procedure based on the following probability 
theory is employed. The procedure will be summarized in the following. 

20 What can be observed by MEG or EEG is several to several hundreds of 

magnetic fields (electric fields) existing near the surface of the brain. In order to 
approximate the current distribution on the virtual curved surface, lattice points are set 
on the virtual curved surface, and a current dipole (or an appropriate current source 
model) is allocated to each lattice point. In order to estimate the current distribution 

25 with high resolution, it is necessary to increase the number of lattice points to increase 
the density of the lattice points. 

When the number of lattice points on the virtual curved surface is increased to be 
larger than the number of observation points, a unique solution cannot be determined. 
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When estimation points is larger in number than the observation points, the number of 
parameters to be estimated becomes larger than the number of observed data, and hence, 
the observed magnetic field would be better recovered even on a virtual surface that is 
positioned deeper than the current source. 
5 In order to cope with such problems, the current distribution on the virtual 

curved surface is estimated utilizing Bayesian estimation theory. At the time of 
Bayesian estimation, prior distribution that represents prior information of the current 
source is used. Specifically, as it is considered that brain current sources exist localized 
at a plurality of positions, a prior distribution that leads to the smallest possible 

10 expansion of current distribution is introduced. Namely, a prior distribution is 

introduced in which a current dipole on a lattice point of which magnitude cannot be 
well identified only from the observed data comes to have a magnitude close to zero. 
This can be realized by introducing hierarchical prior distribution referred to as 
"Automatic Relevance Determination" prior distribution (Reference: R. M. Neal, 

15 Bayesian Learning for Neural Networks, Springer Verlag, 1996). This prior 

distribution will be hereinafter referred to as "localized condition prior distribution." 
The manner how to select a prior distribution introducing prior information other than 
the localized condition will be described later with reference to the second embodiment. 
It is impossible, however, to analytically calculate posterior probability 

20 distribution from the localized condition prior distribution and observed data. 

Therefore, in the present invention, variational Bayes method is used as will be described 
later (Reference: H. Attias, Proc. 15th Conference on Uncertainty in Artificial 
Intelligence pp. 21-30, 1999 and Masa-aki Sato, Neural Computation, 13, 1649-1681, 
2001). It is noted that other method of approximation such as Monte Carlo method 

25 may be used. 

By Bayesian estimation using localized condition prior distribution, it becomes 
possible to obtain a current distribution on the virtual curved surface that recovers the 
observed data and has smallest possible expansion. Further, by comparing model 
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posterior probabilities when estimations are made using virtual curved surfaces of 
different depths, the depth of the current source can be estimated. 

The logarithm of the model posterior probability can be represented as a sum of 
log-likelihood term and model complexity term. The log-likelihood term becomes 
5 larger as recovery error becomes smaller. 

The model complexity term becomes larger when the number of effective current 
dipoles (that is, having a magnitude not smaller than an appropriate threshold) on the 
lattice points becomes smaller. As already described, the recovery error and the 
expansion of the current distribution on the virtual curved surface (number of effective 
10 current dipoles) become the smallest when the depth of the virtual curved surface 
matches the current source. Thus, it can be understood that the model posterior 
probability becomes the largest at this time. In other words, the current source exists 
on the virtual curved surface at a depth at which the model posterior probability 
becomes the largest. 

15 In summary, the present invention provides a method of estimating the position 

of the brain current source, depth direction inclusive, from the observed data of MEG 
and EEG. 

The basic idea of the present invention comes from the fact that an 
electromagnetic field generated by a current source can be recovered by causing an 

20 appropriate current flow over a curved surface existing between the current source and 
an observing surface, and that the expansion of current distribution on the curved 
surface becomes smaller as the curved surface comes closer to the current source, as 
described above. The present invention is characterized in that, base on this idea, the 
position of the current source including the depth direction is estimated by considering 

25 the fact that the model posterior probability becomes the largest when the curved 
surface encompasses the current source in Bayesian estimation of the current 
distribution on the curved surface that recovers the observed data, that is, by considering 
the model posterior probability. 
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(Where there are plurality of current sources) 

Though description has been made on one current source, the method is also 
applicable even when there are a plurality of current sources. 

The electromagnetic field generated by a current attenuates in inverse proportion 
5 to the square of distance, and therefore, the current source closest to the brain surface 

has the largest influence on the observed magnetic field on the brain surface. Therefore, 
it is possible to identify the current sources one by one in order, starting from the one 
closest to the brain surface. 

When the virtual curved surface is moved deeper from the brain surface, the 
10 model posterior probability attains the relative maximum near a current source closest to 
the brain surface (which will be referred to as the first current source). When there are 
two or more current sources, there will be a plurality of local sets of current dipoles 
corresponding in number to the current sources, in the current distribution on the virtual 
curved surface. 

15 The local set of current dipoles is referred to as localized current distribution. 

A local surface including individual localized current distribution is separated and moved 
in the depth direction to find the model posterior probability. When the local surface 
that corresponds to the first current source is moved, the model posterior probability 
attains to the relative maximum at the depth of the first current source. When other 

20 local surfaces are moved, the model posterior probability never attains to the relative 

maximum at the depth of the first current. Thus, the position of the first current source, 
that is, the position of the current source closest to the brain surface can be identified. 

In order to find the second deepest current source, the local surface 
corresponding to the first current source is fixed, the depth of the remaining local 

25 surfaces are aligned and gradually made deeper. Then the model posterior probability 
attains the relative maximum at the second deepest current source. When the 
individual local surface is moved in the depth direction again, the model posterior 
probability attains to the relative maximum at the depth of the second current source, 
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only when the local surface corresponding to the second current source is moved. In 
this manner, the position of the second current source can be identified. The third and 
the following current sources can also be identified in the similar manner. 

The method is advantageous over the method in which the depth of individual 
5 local surface is moved independently, in that the time for computation can significantly 
be reduced. 

(Method in which resolution is increased gradually) 

According to the method of estimating brain current source of the present 
invention, it is possible to identify position of each of the current sources starting from 
10 the one closest to the brain surface in the above described manner and, in addition, it is 
possible to gradually increase the resolution of current source estimation. 

First, a position of a current source is roughly estimated with a low resolution 
(with a small number of lattice points on the virtual curved surface and a small number 
of sample points in the depth direction). In this stage, the position of a local surface 
15 corresponding to each current source is approximately determined. 

Next, estimation is made with higher resolution. At this time, the area of the 
local surface has become smaller as compared with the original virtual curved surface, 
and hence the resolution is naturally higher when the same number of lattice points are 
used. Accuracy in estimating the current source position can be improved by 
20 increasing the resolution in the depth direction as well. If the current distribution is 

further localized when the resolution is made higher, it is possible to estimate again with 
local surface made smaller and the resolution made still higher. 

On the contrary, if the expansion of current distribution does not much vary even 
when the resolution is improved, it means that the current source expands wide. In this 
25 manner, the resolution can be adjusted in accordance with how the current source 
expands. 

[Specific procedure of current source estimation] 

In the following, specific procedure for identifying the position of a "brain 
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current source" will be described in detail, in accordance with the outline above. 
[(I) Preparation for current source estimation] 

First, preparation for the process procedure for estimating the current source will 
be described. 

5 Specifically, for the estimation of current source, the following procedures must 

be taken in advance. 

(1-1) Determination of the shape of virtual curved surface and current model 
(1-2) Further, in order to estimate the virtual current distribution while moving 

the virtual curved surface, it is necessary to determine sample points on the virtual 
10 curved surface and sample points in the depth direction. 

(1-3) Further, as will be described in detail later, it is necessary to determine in 

advance meta parameters to designate the distribution shape of hierarchical prior 

distribution for estimating current distribution as initial values. 

In the following, the procedure of (1-1) Determination of the shape of virtual 
15 curved surface and current model will be described in grater detail. 

(1-1-1) Determination of the shape of virtual curved surface 

The simplest shape of the virtual curved surface is obtained by regarding the 

brain as a hemisphere and assuming various hemispheres of different radii to be the 

virtual curved surfaces. 

20 When a location where existence of a brain current is highly likely such as the 

cerebral cortex has been known in advance by other measuring method such as 
Magnetic Resonance Imaging: MRI, the shape of the virtual curved surface may be 
determined based on such information. 

Particularly when the shape of the cerebral cortex is known with high precision 

25 and it is considered that the brain current source is non-existent in other regions, the 
current distribution estimation of cerebral cortex points may be performed while 
omitting the search in the depth direction, which will be described later. Further, it is 
also possible to perform the search in the depth direction only in a limited region. 
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In this case, the shapes of virtual curved surfaces having different depths may 
generally differ. It is necessary, however, to determine the shapes not to intersect with 
each other. 

(1-1-2) Determination of current model 

As a current model on the virtual curved surface, let us consider a current dipole 
model. It is also possible to consider other current models. 

Consider appropriate lattice points (sample points) {X n jn=l,..., N} on the virtual 
curved surface. Assume a current dipole of which current intensity is j n on each lattice 
point. Here, the magnetic field formed by the current dipole j n at a lattice point X n on 
an observation point Y\ (i=l, I) on the brain surface is given by the following Biot- 
S avail's expression. 

j„x(Y,-X n ) 

Here, jjl represents magnetic permeability, and by way of example, for a vector 
X n , the expression |X n | represents the absolute value of vector X n . As more accurate 
expression, Sarvas's expression with the brain regarded as a sphere filled with 
conductive solution (Reference: J. Sarvas, Phys. Med. Biol. 32, 11-22, 1987) may be 
used, or numerical solution such as given by finite element method or boundary element 
method may be used, considering detailed structure of the brain. In the following, 
Biot-Savart's expression above will be used, for simplicity of description and 
representation. 

Here, the magnetic field formed by the current dipole {j n |n=l, . . N} on an 
observation point Yj (i=l, I) is represented by the following equation. 

w j„x(Y,-X.) 



-21 - 



Assuming that the direction of the magnetic field observed at the observation 
point Yi is a vector Si, c (c=l, C), component B i9C in the direction of vector Si, c of the 
magnetic field at this point can be given as 



„ ■£ , (i.x(Y,-X.))-S t . 



n=l 



Further, when a local gradient of a magnetic field is to be measured as a 
magnetic field to be observed, a differentiation of the equation above by Yi will be 
observed. 

When the direction of a current dipole at a lattice point X n (position vector) is 
restricted, the current dipole can be represented by the following equation, with a 
possible independent direction of the current dipole being b n ,d(d=l, D). In this 
equation, a case where D = 3 corresponds to a case where there is no restriction on the 
direction. 



jn =2-kd ' b 



n,d 



d=l 



In summary, a magnetic field formed by the current dipole at a lattice point {X n 
|n=l, N} of the virtual curved surface on the observation point Yj can be given as 

N D 
n=l d=l 



G(i,c; n,d) = n 



(b M x(Y,-X B ))-S^ 
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Here, j n ,d is an independent component of the current dipole at the lattice point 
X*. Further, the function G (i, c; n, d) represents the component in Si, c direction of the 
magnetic field formed by the current dipole j n , d at the lattice point X n . 

The problem of estimating current distribution may be considered as a problem 
of estimating a current distribution on the virtual curved surface {j n , d | n=l, N; 
d=l, D} from the observed magnetic field {Bi, c |i=l, c=l, C}. 

The measured value of the magnetic field at the observation point described 
above may be given by the following matrix expression. Here, G is referred to as a 
lead field matrix. When Sarvas's expression or a numerical solution such as obtained 
by the boundary element method is used in place of Biot-Savart's expression, the lead 
field matrix will have different values, and other portions of the following description are 
similarly applicable. 

B = G J 

B = ( B i.c |i = l>"--,I;c=l,---,C):(IxC) dimensional vector 

J = ( J M | n = 1 5 • • • , N; d = 1, • • • , D) : (N x D) dimensional vector 

G = (G(i,c;n,d)| i = 1, .-,1; c = 1,- -.,C; n = 1.-.N; d = 1,. -,D) 
: (I x C) x (N x D) dimensional matrix 

(Current source probability model) 

With the current model determined in this manner, the following "current source 
probability model" is introduced for such current distribution estimation. 
(1-1-3) Current source probability model 

A probability model for the current model on the virtual curved surface 
described above will be considered. 

It is assumed that the observed magnetic field is represented as a sum of the 
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magnetic field formed by the current distribution J on the virtual curved surface and the 
observation noise. Further, it is assumed that the observation noise is Gaussian noise 
having an independent variance a 2 at each measuring point. 

More generally, it may be possible to consider Gaussian noise having a multi- 
5 dimensional normal distribution, in which correlation between noises at respective 

measuring points is represented in the form of a covariance matrix. For simplicity of 
description, an isotropic homoscedastic noise model will be employed in the following. 
Specifically, the observed magnetic field is considered as 

B = G J + £ 

J = J virtual curved surface = ( jn,d I n = 1, • • • , N; d = 1,"',D) 

G = (G(i,c;n,d) |i = c = 1,.»,C; n = 1,--,N; d = 1,-,D) 

: Gaussian noise with each component having independent variance a 2 

The probability distribution for the current model may be considered as follows. 
First, a virtual surface at a specific depth, or a set of local surfaces with the depth 
of each local surface identified, will be represented by M. When a current distribution J 
on the virtual curved surface M is given, the probability P (B | J, 0, M) that the 
observed magnetic field is B is represented by the following equation, where 0= 1/a 2 . 

P(B | J,p,M) = exp[-ip|B-G. J| 2 + -(IC)log(0/27t)] 
p = l/a 2 

25 (1-1-4) Hierarchical prior distribution 

As already described, localized condition hierarchical prior distribution will be 
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used as the prior distribution for the current distribution J on the virtual curved surface 

M. 

The prior distribution for the current distribution J before observation 
(probability of J being realized) is represented by the following equation, under the 
assumption of localized condition hierarchical prior distribution. 

P 0 ( J |o,0,M) = exp[-Ipf> n (fX d ) 2 + Tl lo g(K / 2*)] 

Z n=l d=l ^ n=l 

P„0|x,M) = r(p|l/T,y po ) 

Here, T (. .) represents gamma distribution, which is defined below. In the 
expression below, T (y 0 ) is a gamma function, which is also defined below. 

r(P|b, Yo ) - ^(Yo3/b) y » T^e^'* 

P r(y 0 ) 

r(Y 0 )-f dtt^e- 

Further, a and t are hyper parameters introduced to model the current 
distribution J and prior distribution for inverse variance of noise J3. Hierarchical prior 
distributions P 0 (a | M) and P 0 (x|M) for a and x are 

P 0 (a|M) = nr(a n |a 0 ,y o0 ) 

n=l 

P 0 (x|M) = r(T|T 0 ,y T0 ) 
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In the equations above, meta parameters bar oto, y a o, Jo and bar y^o determining 
the distribution shape of the hierarchical prior distributions P 0 (oc|M) and P 0 (x|M) for 
aand x will be described in detail later. Here, "bar" preceding a name of a variable 
means that the variable has the sign " — " thereabove. 

(1-1-5) Bayesian estimation 

When data B of a magnetic field is observed, the posterior probability 
distribution P (J|B, M) that the current distribution is J can be calculated in the following 
manner, using Bayesian theorem. 

P(J|B,M) = J dp da dx P(J, p,a,x |B,M) 

Here, the posterior probability distribution P (J, p, a, x|B, M) is given as 

P(J, p,a,x|B,M) = 

P(B[J, p,M)P 0 (J|a r p,M)P 0 (P|x,M)P 0 (a|M)P 0 (xlM) 

P(B|M) 

P(B|M) = 

JdJdpdadxP(B|J, p,M)P 0 (J|a,p,M)P 0 (p|x,M)P 0 (a|M)P 0 (x|M) 
Using this posterior probability distribution, an expected value E[J|B, M] is given 

as 

J = E[ J |B,M] = jdJd/3dadrf( J, /3,a,r \ J,M) J 
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Further, P (B|M) represents marginal likelihood of the virtual curved surface 
model M. When the depth of a current source is estimated, a number of models having 
different depths are compared. At this time, it is assumed that there is no prior 
information about the depth. Specifically, in the following, it is assumed that P(M) = 
constant. 

When the observation data B is given, the probability that the model M is a true 
model, that is, the model posterior probability P(M|B), is in proportion to the model 
marginal likelihood P(B|M) under the assumption described above, and hence, the 
following relation holds. 

P(M|B)ocP(B|M) 

(1-1-6) Variational Bayes method 

It is generally impossible to analytically calculate the model marginal likelihood 
when the probability model and the hierarchical prior distribution are given. 

Therefore, as a method of calculating by approximation the model marginal 
likelihood, variational Bayes method is used. It is possible to use other method of 
approximation, such as Monte Carlo method and Laplacian approximation. 

The "variational Bayes method" will be briefly outlined, and specific procedure 
will be described later. 

In order to calculate a true posterior distribution P (J, P, a, t|B, M) by 
approximation, a trial posterior distribution Q(J, P, a, x) is considered. 

Closeness between the two probability distributions P (J, 0, a, x|B, M) and Q(J, 
P, a, t) can be calculated by using K-L distance given by the following expressions. 

KL(Q||P) 

= f dJdpdexdxQ( J, p,a,T)log[ Q( J ' P ' a?T) — ] 
j h vv , , s L P (j,p,a,x|B,M) J 
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= logP(B|M)-F(Q)>0 

F (Q) = J d Jdpdadx Q ( J 5 p, a, x ) x 

P(B [J, P,M)P Q ( J [a, p,M)P 0 (P lx,M)P 0 (a |M)P 0 (x |M) 

Q(J,P,a,x) J 

The K-L distance attains to zero only when the two distributions are equal to 
each other, and otherwise it always assumes a positive value. 

When a free energy F(Q) for the trial posterior distribution Q is defined in the 
expression above, the following inequality results. 

F(Q)<logP(B|M) 

Specifically, the distribution Q(J, P, a, x) that maximizes the free energy F(Q) 
becomes equal to the true posterior distribution P (J, P, a, x|B, M), and at this time, the 
free energy is equal to the marginal log-likelihood. 

According to variational Bayes method, the form of function Q(J, P, a, x) is 
restricted to the following form, to maximize the free energy. 

Q(J,P,a,T)=Q J (J,P)Q a (a,T) 

By alternately repeating the step of fixing the second term Qoc on the right side 
of the equation above and maximizing F(Q) with respect to the first term Qj on the right 
side and the step of fixing the first term Qj and maximizing F(Q) with respect to the 
second term Qoc, a distribution Q* is obtained that attains the relative maximum of free 
energy F(Q). 

[(H) Procedure of current source estimation] 

A specific procedure for estimating the current source after the preparation 
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above will be described with reference to the figures. 

Fig. 4 is a flow chart representing an overall procedure of the brain current 
source estimating method in accordance with the present invention. 

Referring to Fig. 4, first, when the process of estimating a position of a brain 
current source starts (step SI 00), values of meta parameters for designating the shape of 
hierarchical prior distribution for estimating the current distribution and an initial value 
of a variable to be estimated are determined (step SI 02). 

Thereafter, initial estimation of the current source is performed, using an initial 
resolution, so as to extract candidates of the current source (step SI 04). 

Then, among the current source candidates extracted in this manner, a position 
of a current source that is closest to the brain surface is estimated (step SI 06). 
Thereafter, depths of other current sources is identified successively (step SI 08). 

After the positions of current sources are identified with the first resolution in 
the above described manner, re-estimation of the positions of the current sources is 
performed with the spatial resolution increased (step SI 10), and the process ends (step 
SI 12). 

Processes of respective steps of Fig. 4 will be described in grater detail in the 
following. 

(II- 1) Initial estimation of current source using initial resolution (extraction of 
current source candidates) 

Fig. 5 is a flow chart representing the process of step SI 04 among the steps 
shown in Fig. 4, that is, the process of initial estimation of the current source using an 
initial resolution. 

Referring to Fig. 5, first, sample points in the depth direction with the initial 
resolution are determined to be {Rk|k=l, K}. Current distribution is estimated for 
the virtual curved surface at each depth Rk. Further, based on the current distribution 
estimated in this manner, posterior probability for each depth Rk is calculated. The 
posterior probability corresponds to the free energy value, which will be described later 
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(step S202). 

For the current distribution at the depth R M at which the posterior probability 
calculated in this manner becomes the highest, a relative maximum point of current 
intensity is found (step S204). 

Further, a local surface that encompasses each relative maximum point is 
determined. Assuming that there are L local surfaces, each local surface will be the 
candidate of localized current source (step S206). 

In the following, the process of step SI 02 of Fig.4 and the process of step S202 
of Fig. 5 that follows, will be described in grater detail. 

(II- 1 - 1 ) Determination of meta parameter values for designating shape of 
hierarchical prior distribution for estimating current distribution and initial values of 
variables to be estimated 

As described above, estimation of a current source is performed using variational 
Bayes method. Therefore, the process of step SI 02 shown in Fig. 4 is performed in the 
following manner, with the application of variational Bayes method. 

By way of example, the meta parameters designating the shape of hierarchical 
prior distribution are determined as follows. 

Ypo = 1 (more generally, 0 < y po ) 

y T0 = 0 (more generally, 0 < y x0 ) 

y a0 = 0.1 (more generally, 0 < y a0 ) 

x 0 = k t • V OT (B), K t = 1 (more generally, k t ~ 1) 

V„(B) = -M; i(B, c -B) 2 

1 i=l c=l 
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— 1 1 c 

1 ' ^ i=l c=l 

a 0 =K a iTr(G r G), k = 10 (more generally, 0 < kJ 

N = D • N 

Further, based on the equations above, each variable is initialized in the 
following manner. 

Y p =^I-C + y po 
y x = Yxo + Ypo 

= «0 
T = T 0 

X G =G' G 

(II- 1-2) Specific process of variational Bayes method for estimating current 
distribution 

(1) Calculation of expected values of parameters J, 0 (J- step process) 
20 Here, expected values of parameters J and 3 are calculated. By this process, 

the free energy F(Q) for Qj is maximized. 

By defining a diagonal matrix A as follows, Qj is derived in accordance with the 
following equations. In the following equations, an expected value of a variable is 
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represented by the variable name with " — " attached thereabove. 

A(n,d;n',d') = S^.a, (n, n' = 1,-,N; d,d' = 1.-.D) 

E = Z G +A 
J = E" 1 G'B 

P = Y P [j\ B - G • J | 2 + i J' A J + Ypo x r 1 
Q,(J,P) = Q J (J|p)Q p (P) 

Q, ( J | P) = exp[-|p(J - J)'Z(J - J) + ±log| Z| + lNlog(p/2:i)] 

Q p (P) = r(pip >Yp ) 

(2) Calculation of expected values of hyper parameters, that is, calculation of 
expected values of parameters a, t (process for a-step) 

Following the J step, expected values of hyper parameters a and x are calculated 
In this step, a process is performed to maximize the free energy F(Q) for Q a . 
The procedure can be represented as 

a„=Ya[Yaoa 0 1 +|p|;j',d+i:(2- 1 )(n,d;n ) d)r 1 (n = l,-,N) 

d=l d=l 

^^YxtY^o'+YpoPr 1 

Q a (a ; T) = r(x|x,Y T )n r ( a n |a n >Ya) 

n=l 

(3) Calculation of free energy 
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Using Qj and Q a calculated through the J step and a step as described above, the 
free energy is calculated in the following manner. 



F = LP + H, + Hp + H a + H T 

LP = -I p | B - G • J | 2 - ilr S-^g + 1(1 C) (< logp > -log27t) 



< logp >= log P + vj/(y p ) - logYp 
dfloe Hy)) 

\(/ (y) = ° \\i : digamma function 

dy 

H, = -l[Tr S" 1 A - logl 2" 1 A I - N] - 1 p J'A J 
2 1 1 2 

H p =Y po [log(xp)-Tp +l] + O(Yp,Yp 0 ) 

* (Y> Yo ) s ( log T(y) - yv|/ (y) + y) - (log r(y 0 ) - y 0 log y 0 + y 0 ) 

+ Y 0 0(Y)-logy) 
H T =Yxo[log(x/x 0 )-(x/f 0 ) + l] 

N 



H a=Z Yao[log(a n /a 0 )-(a n /a 0 ) + l] 



n=l 



In this manner, the process from the J step to calculation of free energy 
described above is repeated until the value of free energy F converges. The value of 
free energy F(Q) after convergence gives an approximation of marginal log-likelihood 
log P (B|M). 

Further, model log-posterior probability differs only by the constant from log P 
(B|M), and therefore, the model having the maximum model posterior probability is the 
same as the model of which free energy is the highest, in accordance with the above 
described approximation. 
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By the above described procedure, posterior probability for the depth Rk can be 
calculated. By performing the processes of steps S204 and S206 of Fig. 5 accordingly, 
it is possible to find candidates of current sources localized to respective local surfaces. 

(II-2) Identification of current source closest to the brain surface 
5 Next, the process of step SI 06 of Fig. 4, that is, the process for identifying the 

current source closest to the brain surface among the candidates of current sources 
found in the manner as described above, will be described. 

Fig. 6 is a flow chart representing a process for identifying a current source 
closest to the brain surface. 
10 First, as the initial value, the value of a variable 1 is set to 1 (step S302). Then, 

the variable 1 is compared with a possible maximum value L of variable 1 (step S304), 
and when the variable 1 is not larger than the maximum value L, the process proceeds to 
step S3 06, and when the variable 1 exceeds the maximum value L, the process proceeds 
to step S324 (step S304). Specifically, the process from step S306 to step S322 is 
15 repeated from 1=1 to 1=L. 

In step S306, depth of a local surface other than the 1th local surface is fixed at 
the depth R M calculated in step S204 of Fig 5. 

The value of variable k is set to 1 (step S3 10). Thereafter, the variable k is 
compared with a possible maximum value K of variable k (step S3 10), and when the 
20 variable k is not larger than the maximum value K, the process proceeds to step S3 12, 
and when the variable k exceeds the maximum value K, the process proceeds to step 
S320. 

In step S3 12, the depth of the 1th local surface is first determined to be Rk. 
Thereafter, current distribution is estimated for the set of L local surfaces (step 

25 S3 14). 

Thereafter, the posterior probability (free energy) for this arrangement of the 
local surfaces is calculated (step S3 16). By incrementing the value of variable k by 
only 1 (step S3 18), the process returns to step S3 10. 
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The process is performed for each of the local surfaces having the depths from 
k=l to k=K, and then, the depth Rm(1) of the first local surface that provides the highest 
posterior probability is calculated (step S320). By incrementing the value of variable l 
only by 1 (step S322), the process returns to step S304. 
5 In this manner, among the depths Rm(1) calculated for each variable 1, one closest 

to the brain surface (shallow) is detected, which is denoted by li (step S3 24). 

Through the steps described above, it follows that the hth local surface 
corresponds to the current source closest to the brain surface, the initial estimate value 
of the depth thereof is calculated as RmOi), and the process terminates (step S326). 
10 (H-3) Identification of depth of current source corresponding to each local 

surface 

Figs. 7 and 8 are flow charts representing a process for identifying depth of a 

current source corresponding to each local surface. 

Referring to Fig. 7, the value of a variable s is set to 1 as an initial value (step 
15 S402). Thereafter, the variable s is compared with a possible maximum value L of 

variable s (step S404), and when the variable s is not larger than the maximum value L, 

the process proceeds to step S406, and when the variable s exceeds the maximum value 

L, the process proceeds to step S434 (step S404). Specifically, the process from step 

S406 to step S432 is repeated from s=l to s=L. 
20 In step S406, when the depths of local surfaces identified so far are represented 

as {Rm(1i) 5 .., Rm(1 s )}, the depths of these local surfaces are fixed at {RmOi), Rm(1 s )}, 

respectively (step S406). 

First, it is assumed that 1 is not equal to any of {U,..., l s }, and that 1 belongs to a 

set of { l, L} (step S408). Then, whether all possible values are considered as the 
25 value of variable 1 or not is determined (step S410), and if all the possible values have 

been considered, the process proceeds to step S430. Otherwise, the process proceeds 

to step S412. 

In step S412, the depth of a local surface that is different from the 1th local 
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surface and not any of {U,..., l s } is fixed at Rm 

The value of variable k is set to 1 (step S414). Thereafter, the variable k is 
compared with a possible maximum value K of variable k (step S416), and when the 
variable k is not larger than the maximum value K, the process proceeds to step S418, 
5 and when the variable k exceeds the maximum value K, the process proceeds to step 
S426. Specifically, the process from step S418 to step S424 is repeated from k=l to 
k=K. 

In step S418, the depth of the first local surface is set to Rk.. 

Thereafter, for the set of L local surfaces, current distribution is estimated (step 

10 s420). 

Further, posterior probability (free energy) for this arrangement of the local 
surfaces is calculated (step 422). 

Referring to Fig. 8, in step S426, after the process to k=K is finished in the 
above described manner, the depth Rm(1) of the 1th local surface attaining the highest 
15 posterior probability is calculated (step S426). 

Then, other variable 1 that is not equal to any of {U, l s } belongs to the set of 
{1, L} and is not yet processed is selected (step S428), and the process returns to 
step S410. 

Through the above described steps, among Rm (1) values corresponding to all the 
20 processed variables 1, the value 1 that is closest to the brain surface is calculated and 
denoted by 1 s+ i (step S430). 

Further, s is replaced by s+1 (step S432), and the process returns to step S404. 
When the process ends for all the possible values s, identificaticm of the depths of 
current sources corresponding to respective local surfaces is finished (step S434). 
25 (II-4) Re-estimation with increased resolution 

Figs. 9 and 10 are flow charts representing the process of step SI 10 shown in 
Fig. 4, that is, the process for re-estimating a position of a current source with higher 
spatial resolution. 
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Referring to Figs. 9 and 10, first, numbers of local surfaces corresponding to the 
current sources estimated using the initial resolution are re-numbered so that the 
surfaces have the numbers 1, 2, 3, L starting from the one closest to the brain surface 
(step S502). Then, corresponding to the expansion of current distribution on each 
local surface, the local surface is made smaller (step S504). 

The depth of the local surface Rm(1) calculated by the process up to step S 108 of 
Fig. 4 is denoted as R M oId (l) (step S506). 

New resolution and search width in the depth direction are respectively 
represented as AR and (K L ■ AR). Further, resolution of lattice points on the local 
surface is also increased (step S508). 

The value of a variable 1 is set to 1 as an initial value (step S510). Thereafter, 
the variable 1 is compared with a possible maximum value L of variable 1 (step S5 12), 
and when the variable 1 is not larger than the maximum value L, the process proceeds to 
step S514, and when the variable 1 exceeds the maximum value L, the process proceeds 
to step S534 (step S512). Specifically, the process from step S514 to step S532 is 
repeated from 1=1 to 1=L. 

In step S5 14, the depth of a local surface 1* other than the 1th local surface is 
fixed to Rm o1<! (l 5 )- 

Further, the process from step S518 to step S526 is performed with k=0, ±1, 

±K L 

First, the value of variable k is set to 0 (step S516). Thereafter, the absolute 
value of variable k is compared with the possible maximum absolute value of K L (step 
S518), and when the absolute value of variable k is not larger than the maximum value 
K L , the process proceeds to step S520, and when the absolute value of variable k 
exceeds the maximum value K L , the process proceeds to step S528 (step S518). 

In step S 520, the depth of the 1th local surface is set to (RM° ld (l)+k ■ AR). 

Thereafter, for the set of L local surfaces, current distribution is estimated using 
the present resolution (step S522). 
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Further, posterior probability (free energy) for this arrangement of the local 
surfaces is calculated (step S524). Then, the value k is set to the next one of {±1, 
±K L }, and the process returns to step S518. 

After the above described process is performed until k=±K L) the value k that 
5 provides the highest posterior probability is calculated in step S528, which value is 
denoted by k M . 

Then, RM old (l) is replaced by R M old (l)+k M ■ AR (step S530). The value of 
variable 1 is incremented by only 1 (step S532), and the process returns to step S512. 

When the above described process has been performed until the value of variable 
10 1 attains to the maximum value L and the resolution has reached the final resolution (step 
S534), the process is terminated (step S538). 

When the resolution is not yet the final resolution (step S534), the resolution of 
the lattice points on the local surface and the resolution in the depth direction are 
increased (step S536), and the process returns to step S510. 
15 By the "method of estimating brain current source" as described above, it 

becomes possible to estimate the position of a brain current source, including the 
position in the depth direction, from observation data of MEG (or EEG). Further, 
such estimation in the depth direction is applicable even when there are a plurality of 
current sources. Still further, the method is applicable when the current sources are 
20 localized as in the case of current dipoles or when the current source has an expansion. 
In addition, by the method, it is possible to estimate how the current source expands. 

By additionally performing the process described with reference to Figs. 9 and 
10 after the estimation with the initial resolution, it becomes possible to successively 
increase resolution for estimating a position. This also means that it is possible to 
25 review with the scope of search restricted based on physiological findings or data 
obtained by other observation method. 
[Result of simulation] 

In the following, simulation results will be described, in which current source 
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positions of an assumed model were estimated in accordance with the method of 
estimating brain current source described above. 

Fig. 1 1 is a top view of a magnetic field distribution in the radial direction 
observed on a surface of a hemisphere, assuming that the human brain is a hemisphere 
5 having the radius of 1 0.0 (arbitrary unit). 

Fig. 1 1 shows a magnetic field distribution on a surface of the hemisphere in 
which a single current source is positioned at a radius of 7.0 from the center, as will be 
described later. In the simulation that will be discussed below, noise having the S/N 
ratio of 0.1 is added to the observation data of the magnetic field. 
10 Fig. 12 represents results of initial estimation of current sources using initial 

resolution. 

Fig. 12 shows the result of calculation of current distribution, in which, as 
described with reference to the process of steps SI 02 to SI 04 of Fig. 4, in order to 
perform initial estimation of the current source using the initial resolution, the depth 
15 (radius) is changed stepwise and on the hemisphere of each depth, the current 

distribution that attains the maximum free energy is calculated in accordance with 
variational Bayes method. 

Fig. 12(a) shows the result where the radius is R=5.0, 12(b) shows the result 
where the radius is R=6.0 and 12(c) shows the result where the radius is R=7.0. 
20 Fig. 13 represents results of initial estimation of current sources using initial 

resolution, where the radius is larger (that is, closer to the surface of the brain). 

Fig. 13(d) shows the result where the radius is R=7.5, 13(e) shows the result 
where the radius is R=8.0 and 13(f) shows the result where the radius is R=9.0. 

Fig. 14 represents magnitude of free energy on each virtual hemisphere, obtained 
25 when current distribution that attains maximum free energy is calculated for virtual 
hemispheres at various depths assumed in the brain. 

Referring to Fig. 14, it can be understood that the maximum point of free energy 
exists between the radii of 7.0 and 8.0, when calculation is made assuming that the 
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hemisphere as a whole is a virtual curved surface. 

As can be seen from Fig.. 13(d), one point of relative maximum exists in the 
current distribution on the virtual hemispherical surface having the radius of R=7.5. 

As the current source is initially estimated using the initial resolution in this 
5 manner, on the virtual curved surface having the maximum free energy, a local surface 
that encompasses the relative maximum point of current distribution is calculated, and 
the process for identifying the current source is performed on the local surface. 

Figs. 15 and 16 represent processes for identifying the current source executed 
further on a local surface including a point of relative maximum, with the resolution of 
10 lattice points on the local surface and the resolution in the depth direction increased. 

Figs. 15(a), 15(b) and 15(c) represent current distribution on local surfaces 
where the radius is R=5.0, R=6.0 and R=7.0, respectively. Similarly, Figs. 16(d), 16(e) 
and 16(f) represent current distribution on local surfaces where the radius is R=7.5, 
R=8.0 and R=9.0, respectively. 
15 Fig. 17 represents the free energy calculated with the depth of local surface 

varied. 

Referring to Fig. 17, the free energy as a result of calculation using the local 
surface has the relative maximum value near the radius of R=7, and from the result, it is 
possible to identify that one current source exists at the position of R=7. 
20 The current distribution of the current source at this time is as shown in Fig. 

15(c), and as already described, it can be understood that the current distribution also 
has the narrowest expansion. 

Fig. 18 is a top view of a magnetic field distribution observed on a surface of the 
brain assumed as a hemisphere, when there are two current sources in the brain. 
25 In the example shown in Fig. 18, it is assumed that the current sources are 

positioned at the radii of R=6.0 and R=8.0. 

Figs. 19 and 20 represent results of initial estimation of current sources using 
initial resolution. 
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Specifically, Figs. 19 and 20 illustrates the process (step SI 04) for performing 
initial estimation of current sources using the initial resolution shown in Fig. 4, when 
there are two current sources in the brain. 

Figs. 19(a), 19(b) and 19(c) represent current distribution on virtual 
5 hemispherical surfaces where the radius is R=5.0, R=6.0 and R=7.0, respectively. 

Figs. 20(d), 20(e) and 20(f) represent current distribution on virtual 
hemispherical surfaces where the radius is R=7.5, R=8.0 and R=9.0, respectively. 

Fig. 21 represents radius dependency of free energy, when current distribution 
that attains maximum free energy is calculated for each virtual hemispherical surface of 
10 the whole hemispherical surfaces. 

As shown in Fig. 21, when the free energy is calculated over the entire 
hemispherical surface, it is understood that the maximum value of free energy exists 
between the radii of R=7 and R=8, that is, near the radius of R-7. 5. 

As can be seen from Fig. 20(d), on the virtual hemispherical surface having the 
15 radius of R=7.5, there are two relative maximum points in the current distribution and 
two local surfaces are determined. 

Next, the result of calculation corresponding to the process of step SI 06 shown 
in Fig. 4 will be described. 

Figs. 22 and 23 illustrates the process for calculating the depth of a 1th local 
20 surface Rm(1) attaining maximum posterior probability, to identify a current source 
closest to the brain surface. 

Therefore, Figs. 22 and 23 represent current distribution when a local surface 
corresponding to one current source is moved while the other local surface is fixed on a 
surface that attains the maximum free energy calculated over the entire spherical surface. 
25 Figs. 22(a), 22(b), 22(c) represent current distribution on local surfaces where 

the radius is R=5.0, R=6.0 and R=7.0, respectively. 

Figs. 23(d), 23(e) and 23(f) represent current distribution on local surfaces 
where the radius is R=7.5, R=8.0 and R=9.0, respectively. 
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Fig. 24 represents local surface position (radius) dependency of free energy 
when a virtual local surface is moved as shown in Figs. 22 and 23. 

As shown in Fig. 24, when a local surface corresponding to one current source is 
moved to the position of R=8 while the other local surface is fixed on a surface that 
5 attains the maximum free energy calculated over the entire spherical surface, the free 
" energy is maximized. 

Therefore, from the result shown in Fig. 24, it can be seen that the position of the 
current source closest to the brain surface is at the radius of R=8.0. 

Thereafter, the depth of the other current source is identified. 
10 Here, the depth of the local surface corresponding to one current source is fixed 

at the radius of R=8 0, and the depth of the local surface corresponding to the other 
current source is moved. 

Figs. 25 and 26 represent current distribution when one local surface is fixed on 
a specific surface and a local surface corresponding to the other current source is moved. 
15 Figs. 25(a), 25(b) and 25(c) represent current distribution on local surfaces 

where the radius is R=5.0, R=5.5 and R=6.0. 

Figs. 26(d), 26(e) and 26(f) represent current distribution on local surfaces 
where the radius is R=6.5, R=7.0 and R=7.5. 

Fig. 27 represents local surface position (radius) dependency of free energy 
20 when a virtual local surface is moved as shown in Figs. 25 and 26. 

As can be seen from Fig. 27, the free energy is maximized at the position of R=6. 
It can be understood that by such a process, even when there are two current 
sources in the brain, it is possible to identify the depth of each of the current sources. 
[Second embodiment] 

25 In the following, the second embodiment of the present invention will be 

described, in which, among the processes implemented by software (or partially by 
hardware) of the "brain current source estimating method," "brain current source 
estimating program," and "brain current source estimating apparatus" in accordance 
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with the first embodiment of the present invention described above, the processes from 
"(I) Preparation for current source estimation" to "(3) Calculation of free energy" are 
modified. 

As a background for describing the modifications in accordance with the second 
embodiment, the concept of the first embodiment will be summarized in the following, 
and of the basic concept, portions that are modified in the second embodiment will be 
described. 

(Concept of virtual curved surface estimation in the first embodiment) 

In the first embodiment, as a current model on the virtual curved surface, current 
dipoles on lattice points are assumed, and a process of solving an inverse problem that 
the current distribution J is estimated from observed magnetic field B on the observation 
surface is performed. Here, such an inverse problem is a so called "ill-posed problem" 
and it is difficult to find an exact solution, as (the number of lattice points at which 
current dipoles exist) is generally larger than (the number of observation points). 

In the first embodiment, in order to solve such a problem, the method of 
"Bayesian estimation" has been used. 

Bayesian estimation utilizes the fact that posterior probability distribution when 
virtual surface M and observed magnetic field B are given can be represented by the 
following equation. 

P (J ]B,M) = P ( B ' J ' M > P ( J ' M > 
P(B|M) 

Here, P(J|B,M) is the probability distribution on the virtual curved surface M 
under the condition that the magnetic field B has been observed, and it represents the 
"posterior probability distribution." P(B|J, M) is "data likelihood," P(J|M) is prior 
distribution and P(B|M) is marginal likelihood. 

The marginal likelihood is given by 
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P(B|M) = j*dJP(B|J,M)P(J |M) 



The posterior probability P(M|B) is given by 
p(M|B) = P(B^P(M) ocp(B|M) 

Assuming that there is no prior information related to the depth, P(M) = 
constant, and therefore, the probability that the model M is the true model when the 
observation data B is given is in proportion to the model marginal likelihood P(M|B) 
under the assumption above, as already described with reference to the first embodiment. 

Specifically, in the first embodiment, when the current model is selected, a model 
that attains the highest posterior probability, or the highest marginal likelihood, has been 
selected. 

For such a process, the marginal log-likelihood log P(B|M) represented by the 
following equation is the object of calculation. 

logP(B|M) =< log(P(B| J,M)) >, -KL(P(J|B,M)||P(J|M)) 

= (expected log-likelihood) - (distance of prior/posterior distribution) 
~ -(recovery error) - (effective degree of freedom of parameters) 

As represented by the equation above, that the marginal log-likelihood is 
maximized corresponds to a concept that the error and the effective degree of freedom 
of parameters are minimized, in other words, that the error and the expansion of current 
distribution are minimized. 

In the first embodiment, the problem of maximizing the marginal log-likelihood 
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described above has been solved by the so called "variational Bayes method" in which 
the trial posterior distribution Q is introduced as an approximation of the posterior 
distribution, and the posterior distribution is determined to maximize the free energy (Q). 
(Concept of prior distribution in the first embodiment) 
5 In the first embodiment, as a prior information of the brain current source, a 

localized condition that the current sources exist at a plurality of positions in the brain is 
used, which is represented in the form of hierarchical prior distribution. Specifically, 
hyper parameters a n (n=l, N) are introduced and the form of the hierarchical prior 
distribution given by the following expressions are assumed. 

10 

P 0 (J|a,p)ccexp[-ip|]a n |J n | 2 ] 

^ n=l 

P 0 ( <*n | a 0 , y a0 ) oc exp [ - y a0 a n / a 0 + ( y a0 - 1) log a n ] 

15 From the expressions above, it can be understood that the hyper parameter a* is 

in proportion to inverse variance (the reciprocal of variance) of current J n . When the 
values of variance of current J n and inverse variance P of noise are known in advance, a 
prior distribution may be used in which the value otn is fixed to the value calculated from 
the values of variance of current and inverse variance of noise. The value of variance 

20 of the current, however, is generally unknown, and therefore, in the first embodiment, a 
hierarchical prior distribution P 0 (otn|bar oco y a o) for the hyper parameter oc, is assumed, 
and Bayesian estimation is performed on a*. 

Here, bar On and Yao are meta parameters for determining the shape of 
hierarchical prior distribution, and common values are assumed for every otn When 

25 Bayesian estimation is performed using the localized condition hierarchical prior 

distribution, Bayesian estimation that recovers the observed data with smallest possible 
number of current dipoles is performed as described above. 
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Minimum norm estimation method, which is often used in current source - 
estimation of MEG (Reference: M. S. Hamalainen and R. J. Ilmoniem, Med. & Biol. 
Eng. & Comput., 32, 35-42, 1994) corresponds, from the view point of Bayesian 
estimation, to a process in which the value 0Cn as the prior distribution is made common 
to all the points, and the value is directly determined to be an appropriate value. 

In the minimum norm estimation, total sum of current intensities at all the points 
is to be minimized, and hence points having small current intensities appear in large 
number in the estimated solution. Therefore, this method is disadvantageous in that it 
is difficult to determine whether such small current distributions represent the true 
current distribution or noise derived from estimation error. 

The localized condition hierarchical prior distribution solves this problem of the 
minimum norm estimation method. 

(Modifications in the second embodiment) 

In the second embodiment, prior distribution considering not only the local 
condition but also continuous condition is assumed. When current source distribution 
is to be estimated with high precision with the spatial resolution in the order of 
millimeters, it is necessary to consider continuity of current distribution. Actually, 
neural cells in the cerebral cortex has a columnar structure having the radius of about 10 
mm, and when a neural activity takes place, it is likely that neural cells in an area having 
the radius of about 20 mm fire simultaneously. Considering such points, in the second 
embodiment, a hierarchical prior distribution is assumed that combines the local 
condition and continuous condition, as will be described in the following. 

In order to introduce the continuous condition, an internal variable Z is newly 
introduced, and the current J is represented as the internal variable Z smoothed through 
a smoothing filter W. Here, though not limiting, a Gauss filter may be used as such a 
smoothing filter W. Now, the prior distribution can be represented as 
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P 0 (J|Z,a ) p)ocexp[-ip|;a n (J n -WZ n ) 2 ] 

n=l 

P 0 (Z|X)ocexp[-Ipf> n |Z n | 2 ] 

^ n=l 

In the prior distribution above, integration of Z can be executed easily, as 
5 represented by 

P 0 (J|a,PA) = JdZP 0 (J|Z,a,P)P 0 (Z|^> 

oc exp [--pJ'CA' 1 + W-A" 1 WT'J] 
2 

10 where A and A represent diagonal matrixes having {cXn|n=l, N} and 

{XJn=l, N} as diagonal components, respectively. 

Specifically, assuming the prior distribution P 0 (J|Z, a, p) Po (Z|X) by introducing 
the internal variable Z is the same as assuming the prior distribution that has a 
covariance matrix P" 1 (A -1 + W • A' 1 • W) for the current J. It is easier, however, to 

15 apply variational Bayes method when the prior distribution is represented by using the 
internal variable Z. Therefore, the prior distribution will be represented by using the 
internal variable Z in the following. Further, the value of hyper parameter is not 
known in advance as in the case of hyper parameter ctn, Bayesian estimation is 
performed onctn and^, assuming the hierarchical prior distribution of the following form. 

20 

P 0 ( <*n | a 0n , y a0n ) oc exp [ - y a0n a n / a 0n + ( y a0n - 1) log a n ] 

P 0 ( I^On , YxOn Y^On 'Kn + ( Y*0n " *) kg 1 
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The hierarchical prior distribution extended in this manner encompasses the 
localized condition prior distribution of the first embodiment. Specifically, when the 
smoothing filter matrix W is identically 0, the hierarchical prior distribution combining 
the continuous condition and the localized condition coincides with the localized 
5 condition prior distribution of the first embodiment. 

It is assumed that meta parameters bar oCno, YaOn, bar ^o n and y^on for determining 
the shape of hierarchical prior distribution may depend more generally on the positions 
of lattice points. When there is only the observation data of MEG (or EEG), there is 
no prior information related to these values. Therefore, as in the first embodiment, 
10 these parameters are substantially determined commonly without any dependence on the 
location of lattice points, as will be described in the following. When observation data 
obtained by other observation means is available, however, the values of meta 
parameters may be determined for each lattice point using such information, as will be 
described with reference to the third embodiment. 
15 (Specific calculation procedure of the second embodiment) 

The calculation procedure described in the following is in most part similar to 
that of the first embodiment. Therefore, in the following, description will be made 
mainly focusing on the modifications from the first embodiment. 

First, as in the first embodiment, when a current distribution J on a virtual curved 
20 surface is given and the observed magnetic field is B, a model probability distribution 
P(B|J, P) is given, using logarithmic representation, by the following equation. 

logP(B| J,p) = -Ip| B - G J | 2 + 1(1 .C)log(P/27c) 

observed magnetic field B: (IXC) dimensional vector 
25 current distribution J: R dimensional vector, R=N X D, N: number of lattice 
points, D: number of components 
lead field matrix: G (IX C) X R dimensional matrix 
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[Hierarchical prior distribution] 

The hierarchical prior distribution in accordance with the second embodiment is 
represented by the following equations. 



20 



logP 0 (J|Z j p,a) = ~P(J-WZyA(J-WZ) + ^Xlog(3a n /27t) 



logP(Z| p,^) = -ipZ'AZ + ^f>g(3^ n /2n) 



(A)(n,d;n' s d') = 5 nn .6 dd .a n 
10 (A)(n,d;n',d') = 8 nn .8 d A l 

(n,n':l,"-,N;d,d' D) 



N 



log P 0 ( a) = X [-Yaon a£ a n + (Y a0n - 1) log a n - $ G (y o0n , a 0 ^ ) ] 

n=l 

logP 0 (X) = f; [-y XOn X- J n ^ n + (y XOn - l)l Qg X n - <D G (Yxon.^o!,) 1 



n=l 



15 logP 0 (p |t) = - Ypo tP + ( Ypo - l)logP - (D G ( Yf30 ,x) 
logP 0 (x) = -y^x + (Y t0 " l)logx - 0 G (y^.To 1 ) 
3> G ( Y, x ) = log T ( y) - y log ( yx ) 



Here, as in the first embodiment, P represents inverse variance of noise. 
[Calculation of trial posterior distribution] 

At this time, it is assumed that the trial posterior distribution is represented by 
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15 



20 



the form of a product of Qj, Q z and Q a , such as Q = Qj (J, p) Q z (Z) Q a (a, X, t). 

Then, the free energy F(Q) is maximized successively for each of Qj, Q z and Q a 
First, in the first step, Q z and Q a are fixed, and F(Q) is maximized for Qj, and the 
following equations result. 

Q J (J,P) = Q J (J|P)Q P (P) 

logQ J (J|P)=-|p(J-J)'E(J-J) + Ilog|E| + l(IC)log(p/27i) 
logQ p (P)= -y p p-'P +(Yp - OlogP - <D G ( Yp , P" 1 ) 

Next, in the second step, Qj and Q a are fixed, and F(Q) is maximized for Q z , 
and the following equation results. 



logQ z (Z)=-ip(Z-Z)'Z z (Z-Z) + Ilog|Z z | + lRlog(p/27 l ) 



Further, in the third step, Qj and Q z are fixed, and F(Q) is maximized for Q a , and 
the following equations result. 

Q a (a,\, T) = Q 0 (a)Q,(X)Q T (x) 



logQ a (a) = £ [-Ya„ a^ 1 a n + (y^ - l)loga n - <D G (Ya„» a. 1 )] 



n=l 



N 



io g Q,(X) = £ [- y ta K l K +(Yx» - Oiogx. -o g (y^, x; 1 )] 



n=l 
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logQ T (t) =- Y t t _1 t + (y T - l)logx - <D G (y T , x" 1 ) 



The specific method of calculating unknown numbers in the equations above will 
be described later. 

As in the first embodiment, a trial posterior distribution that attains the relative 
maximum of free energy F(Q) is found again through repetitive calculation in the second 
embodiment. 

Specifically, by the procedure described below, the trial posterior distribution Q 
is calculated through the first to third steps, and this procedure is repeated until the 
value of free energy F(Q) converges. How to determine the constant term at this time 
will be described later. 

(1) Procedure of the first step (referred to as J- step) 

Expected values of current J and inverse variance p, and variable y p are calculated 
in accordance with the following equations. 



2 = G'G + A = 2 G + A 



J = r 1 (G'B + AWZ) 

Yp=(I.C)/2 + R/2 + y po 




B-G-j| 2 +~(J-WZ) / A(J-WZ)-hiz'AZ + — p-'+Y^T 
2 2 2 



(2) Procedure of the second step (Z-step) 

Thereafter, expected value of internal variable Z is calculated as follows. 



E z = A + W'AW 

± z =I + (A-* W'A)W 



(3) Procedure of the third step (a-step) 

Thereafter, the following calculations are made for the hyper parameters ct„, y^, 
_ D 

Yan ~ 2 ^aOn 

=^3| J„ -(WZ) n f+i|; [(E- 1 )(n.d;n.d) + (WE- 1 W) (n,d;n,d)] + y a0ll c& 
^ D 

Y*n ~" 2 + ^ 0n 

= j P| Z» f + \t (2z )(n.d; n.d) + Ywta 

Z Z d=l 

Yt = Ypo .+ Yto 

Y^" 1 =YpoP + YxoV 

The trial posterior distribution Q can be calculated through the procedures (1) to 
(3) above. 

[Calculation of free energy] 

The free energy of the second embodiment is calculated in accordance with the 
following equation. 
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F = LP + H, + H z + H a + H x + Hp + H T 



Expected log-likelihood LP and model complexity terms Hj, H z , Ha, H^, H p , H t 
are calculated as follows. 

1 7TI~ ^ 1„ 1 



LP = --3 1 B - G J | Tr (E 1 G'G) + — (I C) ( < logP > -log27i) 

2 2 2 

<logP>=logp +vj;(Y p )-logYp 

H ; = ]- [log | A S- 1 1 - Tr (A S" 1 ) + R] - 1 p ( J - WZ )' A (J - WZ ) - -!- TrE' 1 W'A W 
2 2 2 

H p = Y P0 [log(P x) - P f + 1 ] + <D(y p , Yp 0 ) 



1 



2 



H^^tloglE- 1 -TrCS^ + Rl-Ipz'AZ 



2 



H « = Z Yao»[log(a n /a 0n )-(a n /a 0n ) + l] 

n=l 

H X = Z Yxo»[log(>: n /X 0n )-(X n /X 0n ) + l] 

n=l 

H T =Y T0 [log(T/x 0 )-(T/T 0 ) + l] 

In the second embodiment also, the posterior distribution P(J|B) corresponds to 
Qj that maximizes F(Q), and the marginal log-likelihood log(P(B)) corresponds to the 
maximized F(Q). 

Through the above described procedure, it is possible to calculate the marginal 
log-likelihood log(P(B)) of a virtual curved surface of a certain depth. 

Other procedures are the same as those of the first embodiment, and therefore, 
description thereof will not be repeated . 

[Meta parameters of hierarchical prior distribution] 

Selection of meta parameters in the calculation of hierarchical prior distribution 
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will be described. 

When observation data obtained by other observation method is not available, 
the values of meta parameters are selected commonly, regardless of the lattice points, as . 
in the first embodiment. 

5 Specifically, bar x 0 , ypo and can be selected in the same manner as in the first 

embodiment. Further, meta parameters bar oto n , YaOn, bar Xo„ and y X o n can be 
represented as 

a 0n = ^On ~ ^0 
XO Y*0n ~ Y\0n = Y a 0 

Here, bar a 0 and y a0 can be selected in the same manner as in the first 
embodiment. 

(Third embodiment) 

15 In the first and second embodiments, the methods of estimating brain current 

sources have been described, assuming that observation data obtained by other 
observation method is not available. 

Recently, however, various methods have been developed to observe neural 
activities in the brain, and it has been made possible to obtain observation data using a 

20 plurality of observing means under the same experimental conditions. Therefore, as the 
third embodiment, modifications from the first and second embodiments will be 
described, where observation data obtained by other observation methods such as MRI, 
fMRI or PET are available. 

First, when MRI observation data for the same person is available, it is possible 

25 to obtain information related to the position of cerebral cortex or positions of other 

portions of the brain from MRI images. In such a case, it is possible to have the shape 
of virtual curved surface assumed in the first and second embodiments conforming to the 
shape of the cerebral cortex or the shape of other portions of the brain. 
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Further, when it is confirmed that a neural activity under certain experimental 
conditions takes place in the cerebral vertex, it is possible to omit searching in the depth 
direction and to estimate current distribution by placing lattice points only on the 
cerebral cortex. 

5 By fMRI or PET, it is possible to obtain information related to the location of 

neural activities. One must be careful, however, in handling such information, because 
what is measured by fMRI is the amount of blood flow and what is measured by PET is 
the amount of radio isotope. Though these are considered to increase at portions 
having higher neural activities, what is measured is not the current intensity derived from 
10 the neural activities. 

It is considered that the degree of activities measured by fMRI and PET are co- 
related to some extent to the current intensity derived from neural activities. It is not 
clear, however, whether the correlation is linear or non-linear. Further, MEG and EEG 
have temporal resolution in the order of milliseconds, while fMRI has temporal 
15 resolution in the order of seconds, and PET in the order of several tens of seconds. 

Specifically, what is obtained by fMRI or PET is an average over a long period of time 
of the measured results of MEG or EEG. 

Dale et al. proposes a method of utilizing the result of observation of MRI, fMRI 
or PET for the current source estimation by MEG or EEG in Reference: A. M. Dale et 
20 al., Neuron. 26, pp. 55-67 (2000). In this method, however, the activity information of 
fMRI or the like is directly applied as variance information of current. 

From the view point of Bayesian estimation in accordance with the present 
invention, this corresponds to direct designation of the value of inverse variance hyper 
parameter otn. As described above, however, the information obtained by fMRI or the 
25 like is not related directly to the current, and from the temporal view point, it is an 
average over a long period of time of the activities measured by MEG or EEG. 

In contrast, in accordance with the present invention, the information obtained 
by fMRI or the like is not directly designated as the variance information of current but 
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designated at the level of meta parameters of hierarchical prior distribution. Meta 
parameters y a0n and y^on are meta parameters that control reliability of information given 
by the hierarchical prior distribution. Specifically, the smaller the valuesy a0n and yxon, 
the lower the reliability, and the higher the valuesy a0n and y^ 0n , the higher the reliability. 

As described with reference to the second embodiment, the hierarchical prior 
distribution in which localized condition and continuous condition are combined is the 
same as assuming P^A" 1 + W ■ A" 1 • W) as the covariance matrix of current of the 
prior distribution, where A and A are diagonal matrixes having {oCn|n=l, N} and 
{^n|n=l, N} as diagonal components, respectively. 

Specifically, when the value ctn is smaller, the current comes to have larger 
variance, and when the value Tin becomes smaller, the covariance component derived 
from the continuous condition of current becomes larger. Meta parameters bar oco n and 
bar X 0n represent expected values of cXn and Xn when there is no observation data 
obtained from MEG (or EEG). 

When it is expected that neural activity is vigorous at a current lattice point n 
from fMRI or other observation method, the values of meta parameters bar oton and bar 
are made smaller, and bias information may be entered to facilitate the current at this 
current lattice point to assume a larger value. 

In this manner, by introducing measured amount or information obtained by 
other observation method having different temporal scale at the level of meta parameters 
of hierarchical prior distribution, it becomes possible to enter information of other 
observation means as ambiguous information. 

(Simulation results) 

Simulation results in accordance with the second and third embodiments will be 
described. 

In the following, an area including a certain visual cortex at a back portion of 
cerebral cortex obtained by MRI image of a person is used as a virtual curved surface, 
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and the search in the depth direction is omitted. 

As for the simulation data, currents having the current intensity of (1.0, 1.0, 0.5) 
(arbitrary unit) and currents having the current intensity of (1.0, 0.0, 0.5) are caused to 
flow to areas VI, V2 and V5 in the visual cortex, magnetic fields formed by these 
5 currents on the MEG sensor are calculated, and Gaussian noise having the S/N ratio of 
0.1 is added thereto, to provide the observation data. 

Figs. 28 and 29 represent results of Bayesian estimation considering localized 
condition only. In the figures, a vertex of each triangle represents a lattice point 
prepared on the virtual curved surface, and a current dipole is allocated to each lattice 
10 point. The interval between each of the lattice points is about 3 mm in average, and 
spatial resolution is high. Further, the virtual curved surface has a complicated three 
dimensional shape as it is a separated part of cerebral cortex. The figures are two- 
dimensional projection of the three dimensional shape. In the figures, reference 
characters VI, V2 and V5 denote the areas VI, V2 and V5 of the visual cortex. In the 
15 figures, black circles represent points having high estimated current intensity. 

Figs. 30 and 31 represent result of Bayesian estimation considering both 
continuous condition and localized condition. The results of estimation recover the 
position and expansion of current distribution of almost true values. As can be seen 
from Figs. 28 and 29, when Bayesian estimation is made using localized condition, the 
20 position of the current is correctly estimated, whereas the expansion of the current is 
estimated to be too small. From this simulation, it is appreciated that introduction of 
continuous condition is effective when spatial resolution is high. 

Next, a simulation was made with the S/N ratio of noise increased to 0.2 and the 
current distribution unchanged. In that case, even when hierarchical prior distribution 
25 combining the localized condition and continuous condition was used, the expansion of 
current distribution could not be recovered, because of the influence of severe noise. 

Thereafter, a simulation was made assuming that information of observation data 
obtained by other observation means is available under the bad condition as above. In 
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the simulation, it was assumed that information indicating that activities in areas VI, V2 
and V5 are vigorous is obtained from the observation data of fMRL Considering the 
fact that the active areas of fMRI do not always correspond to current activities and that 
fMRI corresponds to an average over a long period of time of MEG observed data, it 
5 was assumed that the areas indicating high activities in fMRI were areas VI, V2 and V5 
expanded by about 20 mm in radius. The same fMRI information was used when the 
current intensities at VI, V2 and V5 were (1.0, 1.0, 0.5) and (1.0, 0.0, 0.5), respectively. 

By the setting above, a situation is simulated in which the active area of fMRI 
does not always correspond to the current intensity, though it is correlated to some 
10 extent. Under such setting, meta parameter of lattice points corresponding to the 
active area of fMRI was set to bar Oo n = bar ?io n = 2 X Tr(G' • G) /tilde N, and meta 
parameter of other lattice points was set to bar a 0n = bar Xo n = 50000 X Tr(G' ■ G)/tilde 
N. Here, the term "tilde" preceding a variable means that the variable has the sign "~" 
thereabove. 

15 The values of y a on and yxon were commonly set toy a o n = Yxon = 0. 1 for all points. 

Figs. 32 and 33 represent the results of estimation at this time. From these 
figures, it can be understood that even when there is considerable noise, correct 
estimation is possible by adding information obtained by other observation means. 
Although the present invention has been described in detail, it is clearly 
20 understood that the same is by way of illustration only and is not to be taken by way of 
limitation, the spirit and scope of the present invention being limited only by the terms of 
the appended claims. 
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